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SUMMARY 


^This  Memorandum  demonstrates  a superposition  method  which  provides  a 
practical  and  efficient  method  for  the  solution  of  linear  elastic  crack  problems 

The  method  approximates  the  linear  elasticity  solution  over  a substructure 
containing  the  crack  tip  by  the  superposition  of  appropriate  singular  stress 
fields  and  a relatively  coarse  finite  element  mesh.  A hybrid  variational 
principle  is  employed  to  ensure  displacement  compatibility  between  this  sub- 
structure and  adjoining  ones  employing  finite  elements  alone. 

A key  feature  of  the  method  is  that  it  may  be  implemented  by  adding 
routines  to  existing  finite  element  packages.  The  particular  implementation 
discussed  here  is  intended  for  use  in  the  aircraft  industry . — - — — 
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I INTRODUCTION 

Linear  theory  predicts  singular  stresses  at  crack  tips1  which  create  serious 

consequences  when  solutions  are  sought  by  numerical  methods  such  as  finite  element 

2 

analysis.  It  is  known  that  the  energy  convergence  of  conventional  finite 
elements  on  refining  the  mesh  is  limited  by  the  order  of  the  singularity,  and  that 
it  is  not  improved  by  the  use  of  higher  order  elements.  As  a consequence, 
although  there  is  still  a limited  amount  of  new  work  employing  conventional 

o i 

analyses'*  ’ , much  of  the  research  effort  in  this  area  has  been  aimed  at  incorpora- 
ting the  singular  stresses  directly  into  the  finite  element  formulation.  This  is 

usually  done  for  a single  element  embedding  the  crack  tip^*^  or  possibly  for  a 

7 8 

group  of  elements  adjacent  to  the  crack  tip  ' . The  resulting  super-elements  are 
capable  of  modelling  the  stresses  at  the  crack  tip  with  considerable  accuracy, 
while  conventional  elements  are  used  to  idealize  the  remaining  structure.  However 
to  achieve  overall  accuracy  further  conditions  must  hold.  Firstly  there  must  be 
adequate  displacement  compatibility  between  the  super-element  and  the  surrounding 
conventional  elements.  Secondly  those  surrounding  elements  must  themselves  be 
small  enough  to  give  a close  representation  of  the  stresses  adjacent  to  the  super- 
element. This  situation  can  lead  to  some  conflict  in  the  choice  of  dimensions  for 
the  super-element.  It  should  be  large  enough  to  cover  the  region  of  large  stress 
gradients  dominated  by  the  1/Vr  term,  yet  small  enough  to  interface  with  a 
sufficiently  fine  exterior  mesh.  If  it  is  necessary  to  refine  the  exterior  mesh 
the  super-element  is  reduced  in  size,  leaving  still  larger  stress  gradients  to  be 
idealized  by  the  conventional  elements.  Thus  the  energy  convergence  using  the 

2 

super-element  is  no  hotter  than  that  achieved  by  the  conventional  finite  elements 
in  the  presence  of  the  singularity. 

This  Memorandum  discusses  a method  which  restores  the  full  rate  of  converg- 
ence to  the  finite  element  scheme  by  enabling  the  singular  functions  to  be  used 
over  a fixed  region  about  the  crack  tip,  while  interfacing  with  an  arbitrarily 
fine  conventional  finite  element  mesh  on  the  exterior.  The  interior  region, 
which  is  treated  here  as  a substructure,  represents  stress  and  displacement  fields 
by  the  superposition  of  finite  elements  and  singular  trial  functions.  The  finite 
elements  are  only  required  to  model  discrepancies  between  the  exact  solution  and 
the  singular  trial  function,  and  are  most  significant  at  the  boundary  at  the 
substructure. 

Although  the  method  gives  correct  convergence  with  mesh  refinement,  it  is 
also  iisportant  that  it  be  competitive  with  super-elements  for  the  coarsest  meshes 
if  it  is  to  be  of  value  in  an  engineering  context. 


■i I ■ 
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The  theoretical  development  of  the  method  is  described  in  Ref  9.  Here  we 
consider  a formulation  of  the  method  which  is  suitable  for  use  with  commercial 
finite  element  packages.  In  particular  the  discussion  relates  to  an  RAE  program 
currently  installed  on  both  the  RAE  and  British  Aerospace  ICL  computers. 

Also  sufficient  numerical  results  are  presented  to  enable  the  technique  to 
be  evaluated  as  a practical  tool  for  handling  singularities  in  linear  elasticity. 
First,  however,  the  method  itself  is  outlined. 

2 METHOD 

The  underlying  idea  of  this  work  is  that  the  accuracy  of  finite  element 
analysis  in  the  presence  of  stress  singularities  may  be  restored  by  the  inclusion 
of  singular  trial  functions  appropriate  to  the  problem.  In  this  instance  the 
functions  correspond  to  stresses  at  a crack  in  a remotely  loaded  infinite  sheet. 

It  would  be  possible  to  include  these  extra  trial  functions  in  sufficient  numbers 
to  n»ake  the  finite  elements  redundant.  Here,  however,  the  policy  is  to  retain 
the  flexibility  of  the  finite  element  method  by  using  the  minimum  numbers  of 
singular  fields,  and  further  restricting  their  use  to  a finite  region  about  the 
crack  tip. 

A variational  principle  is  used  to  minimise  the  difference  between  the  exact 
solution  and  the  approximation.  In  this  instance  the  superposition  of  a singular 
trial  function  and  a mesh  of  simple  displacement  finite  elements  are  used  to  form 
the  approximation.  There  is  some  difficulty  in  maintaining  displacement  compati- 
bility between  these  augmented  trial  functions  and  the  regular  finite  elements  of 
adjacent  substructures.  This  is  overcome  by  using  a hybrid  variational 
principle1®,  and  seeking  the  stationary  point  of  the  functional 

nc  - - IHo1)  + J T^ds  - j TQds  , (1) 

3V  CT 

where  U(o*)  is  the  total  strain  energy  summed  over  all  elements,  o*  and  T* 

are  respectively  stress  fields  and  the  corresponding  tractions  used  as  trial 

functions  in  the  interior  of  each  element,  0 are  displacement  trial  functions 

defined  on  the  boundary  3V  of  each  element,  and  T are  tractions  prescribed 

on  C_  . 

T 

The  finite  element  used  is  the  constant  strain  triangle  whose  piecewise 

p 

linear  displacements  are  denoted  u . The  continuous  functions  used  to  augment 
the  solution  are  denoted  u*  where  the  suffix  i identifies  the  particular 
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function.  In  order  to  facilitate  the  enforcement  of  displacement  continuity 
across  the  interface  Cg  between  the  local  substructure  and  adjacent  ones  it  is 
also  convenient  to  introduce  the  finite  element  field  ut  which  takes  the  same 
values  as  the  corresponding  augmenting  field  ut  at  the  element  nodes.  The 
displacements  defined  on  the  perimeter  of  each  element  are  then  taken  to  be  of 
the  form 


f uF  + £ cu(u*  “ in  t*ie  special  region  and  on  C ^ 


u 


on  Cc 


(prescribed  on  CR) 


(2) 


The  unknowns  represent  the  nodal  connection  quantities  and  the  coefficients 
associated  with  each  augmenting  trial  function.  This  choice  ensures  that  the 
additional  trial  functions  take  zero  values  at  every  node,  enabling  the  use  of 

p 

nodal  connection  quantities  q for  enforcing  compatibility  between  adjacent 
substructures  in  the  usual  manner. 

The  interior  stress  fields  used  are  constructed  by  summing  the  constant 
stresses  over  each  element  and  the  singular  fields  oua*  corresponding  to  the 
displacement  fields  a.u^  . Variation  of  (1)  with  respect  to  the  trial  stress 
fields  yields  the  exact  elasticity  solution  corresponding  to  the  prescribed 

boundary  displacements  whenever  the  element  perimeter  does  not  coincide  with  Cc 

F + ** 

or  Cv  . Thus  the  constant  fields  o and  oT  are  introduced  explicitly  and 

* F 1 + 

correspond  to  the  boundary  displacements  u and  uT  respectively.  The  form 

used  for  the  interior  stress  fields  is  then  given  by 


I 

o - 


oF  + 


l a.  (a*  - at  ♦ a?) 
(•  l'  1 x 1' 


(3) 


The  correction  terms 


oV  are  required  only  for  the  elements  adjacent  to  the 


interface  C_  or  kinematic  boundaries 


Cv  , where  the  augmenting  stress  fields 

Iv 


do  not  correspond  to  the  piecewise  linear  displacements  enforced  on  the  boundary. 
Substitution  of  terms  (2)  (3)  into  the  functional  (1)  gives 


I 


6 


u(uF)  - f u 


FTds 


l ai  - 2U(uF»ui)  + l °jU(ui'uj)  + / ds  - J (UI  - ui) 

^ i ' J ' ' r p ' ' 


Tds 


1 ? °j  /(“*- 2ui) Tjd*  * / uFT*ds'  / (u* -«!)*• 

Jr'  r if  r j.r  ' 


CS+CK 


VCK 


+ \ a. 

j J 


-i 


/ ',iTjd*  * j (“‘  ' Ui)Tjd*  ' U(ai’°j  ) 


CS*CK 


V°K 


- J (u*  " “i  Tj“S) 


CS+CK 


(4a-m) 


Values  for  the  correction  terms  oT  may  be  derived  algebraically  on  an 

1 Q 

element  by  element  basis.  Variation  of  (4)  with  respect  to  ck  gives 


C 
o . 
m 


1 E 


A ..  2 

(1  - v 


- / (u*.  - ut.) 

) J \ nl  nl/ 


ds 


c 

1 Ev 

/ (“»i 

1 

- “Ii)d* 

0 • 
SI 

A n 

(1  - v ) 

c 

1 E 

/(•:. 

- u+.Jds 
•l/ 

nsi 

A 2(1  ♦ v) 

(5) 


J 


where  uY  “ 


(dni-V-,n.i) 


and 


A - l(i,  ♦ l2)h  is  the  area  of  the  element  shown  in  Fig  I.  If  two 
sides  of  the  element  coincide  with  the  boundary  Cg  + contributions  from  each 
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c 

side  are  summed  to  give  o.  . Having  determined  these  unknowns,  the  functional 

1 . F 

(4)  may  be  reduced  to  a quadratic  polynomial  in  cu  and  q . Minimisation  may 
then  be  achieved  by  the  solution  of  a set  of  linear  equations  in  the  normal 
finite  element  manner. 

3 IMPLEMENTATION 

3. 1 General  framework 

The  implementation  described  generates  a suitable  substructure  for  inclusion 
in  commercial  finite  element  packages.  Its  stiffness  matrix  is  presented  in  the 

, form  the 

usual  stiffness  matrix  assembled  from  constant  strain  elements,  while  P , P 

I D 

are  the  corresponding  load  vectors.  The  main  package  would  then  condense  out  the 
freedoms  specified  as  being  internal  giving  the  usual  reduced  matrices 

*BB  " KBB  K3IK1IKIB 

PB  ‘ PB-*B*1PI  • <6> 

The  present  method  however  supplies  extra  rows  and  columns  ^Aj»  kAb*  KAA 
stiffness  matrix  and  terms  P^  in  the  load  vector  corresponding  to  the  additional 
unknowns  ou  . For  computational  convenience  these  unknowns  are  treated  as  the 
freedoms  associated  with  an  extra  node.  Although  these  last  freedoms  will  be 
local  to  the  particular  substructure  it  may  nonetheless  be  convenient  to  treat 
them  as  boundary  freedoms  associated  with  the  highest  level  of  substructuring. 

In  this  way  physically  meaningful  quantities  such  as  the  stress  intensity  factors 
may  be  determined  without  further  processing  of  constituent  substructures. 

The  following  subsections  show  how  the  terms  of  functional  (4)  may  be 
assembled  into  the  augmented  'stiffness'  matrix. 

3.2  Contour  integrals 

Considering  the  terms  of  functional  (4),  the  first  two  terms  only  are  found 
in  conventional  finite  element  formulations,  whilst  several  of  the  extra  terms 
require  integration  over  the  boundary  of  the  substructure. 

Within  the  RAE  program,  the  necessary  values  of  u*  and  T*  for  the 
augmenting  trial  functions  are  supplied  by  a subroutine,  and  the  integrations 
performed  by  six-point  Gauss  quadrature.  The  program  accepts  data  corresponding 
to  each  segment  of  the  boundary,  defined  in  terms  of  its  end  nodes.  The  type  of 


partitioned  form  shown  in  Fig  2.  The  submatrices  Kjj,  K^,  Kj.^ 


boundary  condition  is  specified  and  the  appropriate  values  for  prescribed 
tractions  T or  displacements  u are  read.  If  no  boundary  conditions  are 
defined  at  the  perimeter  of  the  substructure,  a default  option  supplies  interface 


boundary  C conditions. 
The  terms 


J uFT*dt 

CT 

f uFT*c 

J 1 


c +c„ 

K S 


contribute  to  the  submatrix  (KAJ,  K^)  . The  symmetric  part  only  of  the  terms 


J (“*  " 2ut  j T*ds 


[ u*T*c 
J i J 


CS+CK 


is  used  to  form  the  square  submatrix  . The  augmenting  part  PA  of  the  load 

vector  is  formed  from 


/ « - “D 


Tds 


Another  integration  performed  at  the  same  point  in  the  program  is 


I (“i  " UP 


VCK 


which  is  required  for  the  evaluation  of  the  stresses  from  equation  (5), 


correction 


o.  are  stored  for  later  use. 
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The  nodal  loads  corresponding  to  these  stresses  in  each  element  also 
contribute  to  the  submatrix  (K^,  K^)  though  the  term 


/ - “t)TFds 


(4i) 


CS+CK 


The  final  integration 


/ 


F- 

u Tds 


(4b) 


to  determine  the  usual  finite  element  load  vector  (P  , P ) is  simplified  by  the 

IB 

assumption  that  T is  piecewise  constant. 


3.3  Area  integrals 

The  simplicity  of  the  elements  reduces  the  area  integrations  to  simple 


simulation  over  the  elements  of  the  substructure.  Expression  of  the  first  term  in 
the  form 


U(uF)  - J qFKqF 


(4a) 


is  a feature  common  to  all  finite  element  programs.  Here  the  unaugmented  stiff- 
ness matrix  K is  assembled  from  element  stiffness  matrices  which  have  been 


calculated  algebraicly  rather  than  by  numerical  quadrature.  The  contributions  to 
[KIA-  'SJ  resulting  from  the  terms  U(uF,  u^)  are  also  evaluated  on  an  element- 
by-element  basis.  The  corresponding  stresses  ot  are  added  to  the  storage 

n 1 


containing  and  summation  of  the  strain  energy 


u(ot  + o?,ot  ♦ o?) 
1 i j j' 


over  the  elements  of  the  substructure  gives  the  contributions  to  arising 


from  the  symmetric  part  of  terms 


U(o^) 


(4*) 


i A 


10 


J (ui  - ui)T5ds  - 2u(°i»05) 


CS+CK 


(4m) 


/ (UJ  ~ ut)Ttds  - 2U (cK.op 
CS+CK 

U(ut,up 


(4k) 


(4d) 


in  (4). 

This  completes  the  assembly  of  the  augmented  substructure  stiffness  matrix. 
Additional  organisational  routines  may  be  employed  to  rotate  the  reference  frame 
from  the  local  two-dimensional  coordinates  and  embed  it  in  the  globally-defined 
three-dimensional  system. 

4 RESULTS 

The  validation  of  the  superposition  method  in  Ref  9 is  performed  using  a 
single  problem;  namely  a centre  cracked  square  sheet  under  uniaxial  end  load. 

The  points  established  are  that  convergence  of  the  stress  intensity  factor  with 
mesh  refinement  appears  to  be  linear,  and  that  greatest  accuracy  is  obtained  by 
employing  the  singular  functions  over  intermediately  sized  regions.  The  method 
is  also  found  to  be  particularly  accurate  when  applied  to  such  a problem.  Two 
alternative  approaches  for  determining  stress  intensity  factors  are  employed. 

The  first  bases  the  estimate  directly  on  the  coefficients  of  the  singular  trial 
function,  while  the  second  uses  the  program  to  provide  energy  estimates  which  are 
then  used  to  give  stress  intensity  factors  by  strain  energy  release.  The  latter 
is  marginally  more  accurate,  but  in  each  case  accuracy  of  better  than  2%  is 
achieved  for  an  absolutely  regular  mesh  in  which  the  element  size  is  also  the 
semi-crack  length. 

The  test  cases  presented  here  are  more  stringent  insofar  as  they  include 
cracks  growing  from  regions  of  high  stress  gradients  or  towards  reinforcement. 

In  the  first  instance  the  finite  element  mesh  employed  is  only  just  adequate  for 
the  analysis  of  the  uncracked  structure;  even  so  reasonably  accurate  estimates 
for  the  stress  intensity  factor  are  obtained. 

These  results  for  symmetric  cracks  growing  from  a circular  hole  are  compared 
with  those  given  by  Rooke  and  Cartwright*1  in  Fig  3.  The  points  plotted  are 
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obtained  by  strain  energy  release,  based  on  analysis  results  using  the  finite 
element  meshes  shown  with  a/b  *>  0.5,  0.6,  0.7,  0.8,  0.9  . Estimates  obtained 
directly  from  the  values  of  the  coefficients  a are  in  good  agreement. 

The  same  finite  element  idealizations  are  then  used  to  determine  the 
effect  of  various  forms  of  loading  applied  to  the  hole,  this  providing  a prelim- 
inary study  of  cracks  growing  from  pin-loaded  lugs.  Results  shown  in  Fig  4 

corresponding  to  radial  loading  over  the  upper  half  of  the  hole  are  found  to  be 

1 2 

in  close  agreement  (within  2%)  with  those  obtained  by  Parker  using  collocation. 
The  same  figure  shows  that  a point  load  gives  a stress  intensity  factor  up  to  30% 
greater,  the  effect  being  most  marked  for  shorter  crack  lengths.  Fig  4 also 
shows  that  for  these  shorter  crack  lengths  the  effect  of  intermediate  arc  of 
contact  between  pin  and  lug  more  nearly  resembles  the  point  loaded  case. 

The  particular  lug  of  interest  is  shown  in  Fig  5.  The  stress  intensity 

13  . 14 

factors  used  by  Kirkby  and  Rooke  for  this  problem  are  obtained  by  compounding 

known  solutions.  Only  the  results  for  symmetric  cracks  are  checked,  but  here 

agreement  is  sufficiently  close  to  enhance  confidence  in  each  technique.  The 

effects  of  rounding  the  end  of  the  lug  and  of  reducing  its  length  are  also 

considered. 

The  final  problem  applies  the  technique  to  a reinforced  panel  containing  a 
series  of  in  line  cracks.  A regular  12  x 12  mesh  gave  the  points  plotted  in 
Fig  6 against  comparison  curves  from  Rooke  and  Cartwright1 1 which  agree  to 
within  2%. 
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Fig  1 Element  geometry  and  boundary  displacements 
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Fig  2 Partitioning  of  substructure  stiffness  matrix 


Fig  3 Kj  for  symnetric  cracks  at  a circular  hole  in  a rectangular  sheet 
under  end  load 
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Kj  for  symmetric  cracks  at  variously  loaded  hole  in  rectangular  sheet, 
reacted  by  end  load 


Fig  5a  Effect  on  Kj  of  varying  pin  load  and  lug  end  geometry 
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